HW

1. Problem 11.1

data <- read_excel("Chapter11_exercises_data.xls")
## New names:
## • `` -> `...8`
## • `` -> `...9`
data$GSF <- as.numeric(data$GSF)
data$GSJ <- as.numeric(data$GSJ)

data <- data %>% 
  dplyr::select(Date, GSF, GSJ) %>% 
  na.omit()
# Split into estimation and prediction samples
train <- data[1:(nrow(data) - 20), ] %>% dplyr::select(GSF, GSJ)
test <- data[(nrow(data) - 19):nrow(data), ] %>% dplyr::select(GSF, GSJ)

summary(train)
##       GSF               GSJ        
##  Min.   :-4.6800   Min.   :-3.232  
##  1st Qu.: 0.2828   1st Qu.: 0.400  
##  Median : 2.2296   Median : 1.635  
##  Mean   : 2.1558   Mean   : 1.892  
##  3rd Qu.: 3.8616   3rd Qu.: 3.311  
##  Max.   :10.5140   Max.   : 9.799
summary(test)
##       GSF               GSJ        
##  Min.   :-8.7830   Min.   :-8.760  
##  1st Qu.:-4.4226   1st Qu.:-3.623  
##  Median :-0.8497   Median :-1.162  
##  Mean   :-1.7085   Mean   :-1.133  
##  3rd Qu.: 0.5783   3rd Qu.: 1.406  
##  Max.   : 6.4845   Max.   : 4.582

GSF and GSJ represent the quarterly house price growth rate calculated based on the original house price indices for San Francisco-Oakland-Fremont and San Jose-Sunnyvale-Santa Clara, respectively. For the train sample, GSF has a median of 2.2296 and mean of 2.1558. The middle 50% of data lies between 0.2828 and 3.8616. For the train sample of GSJ, the median is 1.635 and mean is 1.892, while the middle 50% of data lies between 0.400 and 3.311. In both test samples, the median and mean are both negative, and the min, max, and quartiles are all much lower compared to the train data. This is because of the dynamics of the time period from 2008 to 2012.

# Create TS objects using estimation sample
GSF_train <- ts(train$GSF,
            start = c(1975, 2), 
            end = c(2007, 3), 
            frequency = 4)

GSJ_train <- ts(train$GSJ, 
              start = c(1975, 2), 
              end = c(2007, 3), 
              frequency = 4)

# VAR model
VAR_data <- window(ts.union(GSF_train, GSJ_train), start = c(1975, 2), end = c(2007, 3))

VARselect(VAR_data, lag.max = 10)
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      3      3      3 
## 
## $criteria
##               1        2        3        4        5        6        7        8
## AIC(n) 1.260763 1.259156 1.003914 1.028296 1.034471 1.071311 1.014798 1.077252
## HQ(n)  1.317364 1.353491 1.135982 1.198098 1.242007 1.316580 1.297802 1.397989
## SC(n)  1.400138 1.491447 1.329121 1.446420 1.545511 1.675267 1.711671 1.867041
## FPE(n) 3.528188 3.522788 2.729665 2.797876 2.816523 2.924191 2.766069 2.947888
##               9       10
## AIC(n) 1.071033 1.080588
## HQ(n)  1.429504 1.476793
## SC(n)  1.953739 2.056210
## FPE(n) 2.934116 2.967940

We will use a lag length of 3 for the VAR model.

var_model <- VAR(y = VAR_data, p = 3)
summary(var_model)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: GSF_train, GSJ_train 
## Deterministic variables: const 
## Sample size: 127 
## Log Likelihood: -418.202 
## Roots of the characteristic polynomial:
## 0.9124 0.6813 0.6813 0.5674 0.319 0.319
## Call:
## VAR(y = VAR_data, p = 3)
## 
## 
## Estimation results for equation GSF_train: 
## ========================================== 
## GSF_train = GSF_train.l1 + GSJ_train.l1 + GSF_train.l2 + GSJ_train.l2 + GSF_train.l3 + GSJ_train.l3 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## GSF_train.l1  0.27804    0.12489   2.226   0.0279 *  
## GSJ_train.l1  0.65073    0.14912   4.364 2.72e-05 ***
## GSF_train.l2 -0.16813    0.13269  -1.267   0.2076    
## GSJ_train.l2 -0.33706    0.17370  -1.940   0.0547 .  
## GSF_train.l3  0.57698    0.12499   4.616 9.88e-06 ***
## GSJ_train.l3 -0.08235    0.15007  -0.549   0.5842    
## const         0.16928    0.21193   0.799   0.4260    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 1.653 on 120 degrees of freedom
## Multiple R-Squared: 0.6244,  Adjusted R-squared: 0.6056 
## F-statistic: 33.25 on 6 and 120 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation GSJ_train: 
## ========================================== 
## GSJ_train = GSF_train.l1 + GSJ_train.l1 + GSF_train.l2 + GSJ_train.l2 + GSF_train.l3 + GSJ_train.l3 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## GSF_train.l1 -0.06316    0.10844  -0.583 0.561319    
## GSJ_train.l1  1.00212    0.12947   7.740 3.48e-12 ***
## GSF_train.l2 -0.19186    0.11521  -1.665 0.098453 .  
## GSJ_train.l2 -0.29834    0.15081  -1.978 0.050200 .  
## GSF_train.l3  0.37086    0.10852   3.417 0.000864 ***
## GSJ_train.l3 -0.01054    0.13030  -0.081 0.935652    
## const         0.28633    0.18401   1.556 0.122327    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 1.435 on 120 degrees of freedom
## Multiple R-Squared: 0.637,   Adjusted R-squared: 0.6189 
## F-statistic:  35.1 on 6 and 120 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           GSF_train GSJ_train
## GSF_train     2.733     1.688
## GSJ_train     1.688     2.060
## 
## Correlation matrix of residuals:
##           GSF_train GSJ_train
## GSF_train    1.0000    0.7112
## GSJ_train    0.7112    1.0000

In the equation for GSF, the second lag of GSF and the third lag of GSJ have p-values higher than 0.05, so we could consider removing them from the model. The equation for GSF has an adjusted R-squared value of 0.6056, which is pretty good. In the equation for GSJ, the first and second lags for GSF and the third lag of GSJ have p-values higher than 0.05, so they may not have a meaningful influence on GSJ. The model has an adjusted R-squared value of 0.6189, similar to the first equation.

2. Problem 11.2

grangertest(GSF_train, GSJ_train, order = 3)
## Granger causality test
## 
## Model 1: GSJ_train ~ Lags(GSJ_train, 1:3) + Lags(GSF_train, 1:3)
## Model 2: GSJ_train ~ Lags(GSJ_train, 1:3)
##   Res.Df Df      F  Pr(>F)   
## 1    120                     
## 2    123 -3 4.0141 0.00921 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(GSJ_train, GSF_train, order = 3)
## Granger causality test
## 
## Model 1: GSF_train ~ Lags(GSF_train, 1:3) + Lags(GSJ_train, 1:3)
## Model 2: GSF_train ~ Lags(GSF_train, 1:3)
##   Res.Df Df      F    Pr(>F)    
## 1    120                        
## 2    123 -3 6.5504 0.0003863 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using an order of 3, the granger-causality tests both have a p-value below 0.05, so they suggest that GSF granger-causes GSJ and that GSJ also granger-causes GSF. This indicates that past values of GSF can be used to predict future values of GSJ, and vice versa.

3. Problem 11.3

plot(irf(var_model))

The first plot shows that a shock in GSF causes an immediate significant response in GSF followed by a gradual decline around the second lag. This suggests that the economy responds significantly to its own shocks in the short term and maintains a positive response throughout. The second plot similarly shows that a shock in GSF causes a positive and significant response in GSJ at first followed by a gradual decline to near baseline level around the second lag, indicating a connection between the two economies. The third plot shows that a shock in GSJ causes a positive response in GSF in the short term followed by a gradual decline to baseline by the third lag, at which point the response is about zero. The final plot shows that a shock in GSJ causes a positive and significant response in GSJ in the short term followed by a gradual decline to baseline by the third lag, at which point the response is roughly zero throughout.

4) Problem 7.8

a)

retaildata <- read_excel("retail.xlsx", skip=1)

myts <- ts(retaildata[,"A3349873A"],
           frequency=12, start=c(1982,4))

autoplot(myts)

Multiplicative seasonality is necessary because the seasonal fluctuations vary with time. The amplitude of this data grows with time, signaling multiplicative tendencies.

b)

fit1 <- hw(myts, seasonal = "multiplicative")
fit2 <- hw(myts, seasonal = "multiplicative", damped = TRUE)

autoplot(myts) +
  autolayer(fit1, series="HW multiplicative forecasts", PI=FALSE) +
  autolayer(fit2, series="HW multiplicative forecasts - damped", PI=FALSE)

c)

rmse1 <- accuracy(fit1)["Training set", "RMSE"]
rmse1
## [1] 13.29378
rmse2 <- accuracy(fit2)["Training set", "RMSE"]
rmse2
## [1] 13.30494

The RMSEs are basically the same, but the one for the multiplicative method is slightly lower so I will prefer that one over the damped method.

d)

tsdisplay(fit1$res)

lbtest <- Box.test(fit1$res, lag=20, type="Ljung-Box")
lbtest
## 
##  Box-Ljung test
## 
## data:  fit1$res
## X-squared = 35.369, df = 20, p-value = 0.01822

No, the residuals from the best method do not exhibit white noise. There are a few lags that cross the bands, signalling autocorrelation. Also, the Ljung-Box test (p = 0.01822) confirms that the residuals are not white noise.

e)

myts.train <- window(myts, end=c(2010,12))
myts.test <- window(myts, start=2011)

fit3 <- hw(myts.train, h=36, seasonal = "multiplicative")
fc <- snaive(myts.train, h=36)

autoplot(myts) +
  autolayer(fit3, series="multiplicative", PI=FALSE) +
  autolayer(fc, series="snaive", PI=FALSE)

accuracy(fit3, myts.test)
##                       ME      RMSE       MAE          MPE      MAPE      MASE
## Training set  0.03021223  9.107356  6.553533  0.001995484  3.293399 0.4107058
## Test set     78.34068365 94.806617 78.340684 19.945024968 19.945025 4.9095618
##                    ACF1 Theil's U
## Training set 0.02752875        NA
## Test set     0.52802701  1.613903
accuracy(fc, myts.test)
##                     ME      RMSE      MAE       MPE      MAPE     MASE
## Training set  7.772973  20.24576 15.95676  4.702754  8.109777 1.000000
## Test set     81.744444 100.00869 82.06667 20.549055 20.669738 5.143067
##                   ACF1 Theil's U
## Training set 0.7385090        NA
## Test set     0.6830879   1.67023
rmse3 <- accuracy(fit3, myts.test)["Test set", "RMSE"]
rmse3
## [1] 94.80662
rmse4 <- accuracy(fc, myts.test)["Test set", "RMSE"]
rmse4
## [1] 100.0087

Yes, the holt-winters multiplicative method beats the seasonal naive approach. The rmse is lower and the graph shows it is slightly closer to the actual test set values than the snaive approach.

5) Problem 7.9

#Box-Cox transformation
lambda <- BoxCox.lambda(myts)
myts.bc <- BoxCox(myts, lambda)
myts.bc <- ts(myts.bc[,1], frequency=12, start=c(1982,4))

#STL decomposition
stl_decomp <- stl(myts.bc, s.window="periodic", robust = TRUE)
autoplot(stl_decomp)

#get seasonally adjusted data
seasonal_component <- stl_decomp$time.series[, "seasonal"]
seasonally_adjusted <- myts.bc - seasonal_component

#define training/test sets
myts_train <- window(myts.bc, end=c(2010,12))
myts_test <- window(myts.bc, start=2011)

#perform ETS
ets.fit <- ets(myts_train)
ets.forecast <- forecast(ets.fit, h = 36)

#check RMSE
rmse5 <- accuracy(ets.forecast, myts_test)["Test set", "RMSE"]
rmse5
## [1] 0.5734974

The RMSE of 0.5734974 is the lowest compared to the other best previous forecasts.

Question 6 11-a)

data(visitors)
autoplot(visitors)

The data exhibits a multiplicative increasing trend. Additionally, a seasonal pattern is evident, recurring every year.

11-b)

fc <- hw(subset(visitors,end=length(visitors)-24),
         damped = TRUE, seasonal="multiplicative", h=24)
autoplot(visitors) +
  autolayer(fc, series="HW multi damped", PI=FALSE)+
  guides(colour=guide_legend(title="Yearly forecasts"))

11-c)

Due to the increase in the magnitude of the visitor numbers, a multiplicative seasonal model is necessary.

11-d)

visitors_train <- subset(visitors, end = length(visitors) - 24)
visitors_test <- subset(visitors, start = length(visitors) - 24)

# 1. ETS model
fc_ets <- forecast(ets(visitors_train), h = 24)
plot(fc_ets)

# 2. Additive ETS on Box-Cox transformed series
lambda <- BoxCox.lambda(visitors_train)
fc_ets_add_BoxCox <- forecast(ets(visitors_train, lambda = lambda, additive.only = TRUE), h = 24)
plot(fc_ets_add_BoxCox)

# 3. Seasonal naive method
fc_snaive <- snaive(visitors_train, h = 24)
plot(fc_snaive)

# 4. STL decomposition on Box-Cox transformed data with ETS on seasonally adjusted data
visitors_train_bc <- BoxCox(visitors_train, lambda)
stl_out <- stl(visitors_train_bc, s.window = 12)  # Use seasonal argument
visitors_sa <- stl_out$time.series
plot(visitors_sa)

11-e)

checkresiduals(fc_ets)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,M)
## Q* = 20.677, df = 24, p-value = 0.6577
## 
## Model df: 0.   Total lags used: 24
checkresiduals(fc_ets_add_BoxCox)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 15.206, df = 24, p-value = 0.9146
## 
## Model df: 0.   Total lags used: 24
checkresiduals(fc_snaive)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 295.02, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The Ljung-Box test reveals that the seasonal naive method fails to capture all the trends in the data. This is because its residuals exhibit a strong autocorrelation, indicated by a very high Q statistic and a very low p-value from the test. On the other hand, the Additive ETS model on Box-Cox transformed series emerges as the best candidate. This model demonstrates the lowest Q statistic and the highest p-value in the Ljung-Box test, suggesting minimal autocorrelation in the residuals and a good fit for the data.

11-f)

# first, make functions to make model to yield forecast class object
fets_add_BoxCox <- function(y, h) {
  forecast(ets(
    y,
    lambda = BoxCox.lambda(y),
    additive.only = TRUE
  ),
  h = h)
}
fstlm <- function(y, h) {
  forecast(stlm(
    y, 
    lambda = BoxCox.lambda(y),
    s.window = frequency(y) + 1,
    robust = TRUE,
    method = "ets"
  ),
  h = h)
}
fets <- function(y, h) {
  forecast(ets(y),
           h = h)
  }

# I'll compare the models using RMSE
sqrt(mean(tsCV(visitors, snaive, h = 1)^2, na.rm = TRUE))
## [1] 32.56941
sqrt(mean(tsCV(visitors, fets_add_BoxCox, h = 1)^2,
          na.rm = TRUE))
## [1] 18.8505
sqrt(mean(tsCV(visitors, fstlm, h = 1)^2,
          na.rm = TRUE))
## [1] 17.49642
sqrt(mean(tsCV(visitors, fets, h = 1)^2, na.rm = TRUE))
## [1] 18.52985
sqrt(mean(tsCV(visitors, hw, h = 1, 
               seasonal = "multiplicative")^2,
          na.rm = TRUE))
## [1] 19.62107

tsCV errors show that the best model is the STL + ETS(M, A, N) model and the worst model is seasonal naive model. If I hadn’t calculated accuracy using test set, I couldn’t have known that the forecasts from seasonal naive method were the most accurate ones.

Question 7 13-a)

data(hsales)
is.ts(hsales)
## [1] TRUE

The data transformation is necessary due to its non-stationary nature (lack of mean reversion).

13-b)

tsdisplay(hsales)

hsales_diff <- diff(hsales, lag = 12 , differences = 1)
tsdisplay(hsales_diff)

Using seasonal difference would take the seasonal pattern present in the data that repeat every 12 months into account.

13-c)

#hsales_diff
hsales_model <- arima(hsales_diff, order = c(13, 0, 0))
AIC(hsales_model)
## [1] 1641.46
hsales_model_2 <- arima(hsales_diff, order = c(25, 0, 0))
AIC(hsales_model_2)
## [1] 1641.232

Based on the exponentially decaying ACF with significant spikes at lags 13 and 25, I’ve chosen to evaluate ARIMA(13,0,0) and ARIMA(25,0,0) models. Since a lower AIC value indicates a better fit, ARIMA(25,0,0) would be a better model.

13-d

coef(hsales_model_2)
##          ar1          ar2          ar3          ar4          ar5          ar6 
##  0.836300050 -0.018559397  0.007202887  0.004160205  0.068567763  0.104486076 
##          ar7          ar8          ar9         ar10         ar11         ar12 
## -0.085498706  0.041765085 -0.079485039  0.093966567 -0.002479969 -0.581887389 
##         ar13         ar14         ar15         ar16         ar17         ar18 
##  0.491495817  0.007722343 -0.014495612  0.024323854  0.061152752  0.006687802 
##         ar19         ar20         ar21         ar22         ar23         ar24 
## -0.072474289 -0.054275164  0.008676560 -0.011662305  0.011350102 -0.259985765 
##         ar25    intercept 
##  0.273731526 -0.189719054
checkresiduals(hsales_model_2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(25,0,0) with non-zero mean
## Q* = 14.475, df = 3, p-value = 0.002325
## 
## Model df: 25.   Total lags used: 28

The residuals exhibit an approximately normal distribution centered around zero, suggesting characteristics close to a white noise process.

13-e)

plot(forecast(auto.arima(hsales),h=24))

13-f)

hsales %>% ets() %>% forecast(h=24) %>% autoplot()

The ETS forecast exhibits a stronger seasonal pattern compared to the other models.

Question 8: 8.18.

a.

gdp.NA <- Quandl("FRED/GDP", api_key= "mVHzLbBisoncyRnYC27C", type = "ts", transform = "rdiff")
str(gdp.NA)
##  Time-Series [1:299] from 1947 to 2022: 0.0115 0.0147 0.0407 0.0231 0.0257 ...
head(gdp.NA)
##            Qtr1       Qtr2       Qtr3       Qtr4
## 1947            0.01153131 0.01470516 0.04070757
## 1948 0.02308803 0.02568281 0.02432063

b. Plot graphs of the data, and try to identify an appropriate ARIMA model.

autoplot(gdp.NA)

ndiffs(gdp.NA)
## [1] 1
ggtsdisplay(diff(gdp.NA))

auto.arima(gdp.NA)
## Series: gdp.NA 
## ARIMA(1,1,3)(2,0,0)[4] 
## 
## Coefficients:
##          ar1      ma1     ma2      ma3     sar1     sar2
##       0.1792  -0.9982  0.2036  -0.1514  -0.0820  -0.1125
## s.e.  0.2982   0.2963  0.2601   0.0632   0.0639   0.0814
## 
## sigma^2 = 0.0001602:  log likelihood = 881.22
## AIC=-1748.44   AICc=-1748.05   BIC=-1722.56

From auto.arima, ARIMA(1,1,3)(2,0,0) model will fit well to the data.

c. Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?

gdp.NA_arima <- auto.arima(gdp.NA)
autoplot(gdp.NA_arima$residuals)

checkresiduals(gdp.NA_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,3)(2,0,0)[4]
## Q* = 6.7797, df = 3, p-value = 0.07926
## 
## Model df: 6.   Total lags used: 9

The residuals are like white noise.

e. Use your chosen ARIMA model to forecast the next four years.

fc_gdp.NA_arima <- forecast(
  gdp.NA_arima, h = 4
)

autoplot(fc_gdp.NA_arima)

The forecast values of gdp production are first increasing. But then increase will stop and decrease a bit, before dampening to constant.

f. Now try to identify an appropriate ETS model.

gdp.NA_ets <- ets(gdp.NA)

gdp.NA_ets
## ETS(A,N,N) 
## 
## Call:
##  ets(y = gdp.NA) 
## 
##   Smoothing parameters:
##     alpha = 0.0586 
## 
##   Initial states:
##     l = 0.017 
## 
##   sigma:  0.0129
## 
##       AIC      AICc       BIC 
## -890.9744 -890.8931 -879.8731

Chosen model is ETS(A, N, N).

g. Do residual diagnostic checking of your ETS model. Are the residuals white noise?

checkresiduals(gdp.NA_ets)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,N)
## Q* = 32.364, df = 8, p-value = 8.014e-05
## 
## Model df: 0.   Total lags used: 8

The residuals are not like white noise.

h. Use your chosen ETS model to forecast the next four years.

fc_gdp.NA_ets <- forecast(
  gdp.NA_ets, h = 4
)

autoplot(fc_gdp.NA_ets)

The forecast values of difference in growth rate in GDP is stationary.

i. Which of the two models do you prefer?

I prefer ARIMA model. Because it has better estimate given the residuals behaving like white noise. Also, I think at the current state, with the US government trying to reduce inflation, there should be slight decrease in the growth rate of GDP.

Question 9

beer=read.csv("beer.csv",header=T,dec=",",sep=";")
beer=ts(beer[,1],start=1956,freq=12) #,1 to take the values
lbeer <-log(beer) #make it has constant variance

train <- window(lbeer, end=c(1988,1))
h <- length(lbeer) - length(train)
ETS <- forecast(ets(train), h=h)
ARIMA <- forecast(auto.arima(train, lambda=0, biasadj=TRUE),
  h=h)
STL <- stlf(train, lambda=0, h=h, biasadj=TRUE)
NNAR <- forecast(nnetar(train), h=h)
TBATS <- forecast(tbats(train, biasadj=TRUE), h=h)
Combination <- (ETS[["mean"]] + ARIMA[["mean"]] +
  STL[["mean"]] + NNAR[["mean"]] + TBATS[["mean"]])/5

autoplot(lbeer) +
  autolayer(ETS, series="ETS", PI=FALSE) +
  autolayer(ARIMA, series="ARIMA", PI=FALSE) +
  autolayer(STL, series="STL", PI=FALSE) +
  autolayer(NNAR, series="NNAR", PI=FALSE) +
  autolayer(TBATS, series="TBATS", PI=FALSE) +
  autolayer(Combination, series="Combination") +
  xlab("Year") + ylab("$ billion") +
  ggtitle("log of Montly production of beer in Autralia")

c(ETS = accuracy(ETS, lbeer)["Test set","RMSE"],
  ARIMA = accuracy(ARIMA, lbeer)["Test set","RMSE"],
  `STL-ETS` = accuracy(STL, lbeer)["Test set","RMSE"],
  NNAR = accuracy(NNAR, lbeer)["Test set","RMSE"],
  TBATS = accuracy(TBATS, lbeer)["Test set","RMSE"],
  Combination = accuracy(Combination, lbeer)["Test set","RMSE"])
##         ETS       ARIMA     STL-ETS        NNAR       TBATS Combination 
##  0.07881405  0.07829543  0.08169119  0.09785364  0.07817191  0.07668909

All of these individual models have really high accuracy, but the combination approach is even better.